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Abstract Three-dimensional numerical simulations have been used to study 
the scattering of a surface-gravity wave packet by vertical magnetic flux tubes, 
with radii from 200 km to 3 Mm, embedded in stratified polytropic atmosphere. 
The scattered wave was found to consist primarily of m = (axisymmetric) 
and m = 1 modes. It was found that the ratio of the amplitude of these two 
modes is strongly dependant on the radius of the flux tube: The kink mode is 
the dominant mode excited in tubes with a small radius while the sausage mode 
is dominant for large tubes. Simulations of this type provide a simple, efficient 
and robust way to start understanding the seismic signature of flux tubes, which 
have recently began to be observed. 

Keywords: Helioseismology, Direct Modeling; Waves, Magnetohydrodynamic; 
Magnetic Fields 



1. Introduction 

Magnetic flux tubes couple the different layers of the solar atmosphere, and 
waves propagating along the magnetic-field lines are a possible source of the 
mechanical-energy transport required to heat the chromosphere and corona. The 
waves also provide a helioseismic means of constraining the properties of the mag- 
netic tubes which is complementary to observatio nal studies with SOHO/MDI 
data ( Duvall Gizon. and Birchil2006f ) and Hinode ( Fuiimura and Tsunetd , l2009[ ).| 



In the latter work the authors found that the observed Poynting flux is sub- 
stantial when compared to the flux required to heat the chromosphere and 
corona. Unfortunately with observations at a single height they were unable to 
identify what mixture of m = 0, axisymmetric sausage-modes, and m = 1, kink- 
modes, they were observing. Theoretical studies of the interaction of flux tubes 
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and waves have bee n performed in the past, fo r example iRoberts and Webb 



(Il979h JWilsonl dl98dh JSpruitl (Il98ll) JSpruitl (Il982h JBogdan and Knolke 



Solankil (Il993h. Bogdan et al. 



(Il996 |), [Hasan and Kalkofenl (Il999h lTi~ 



Ikel 
irrvl 



1989), 



2QQQh. 



Gizon, Hanasoge, and Birchl ( 2006 ) J Jain and Gordovskvvl (|2008f ) jHanasoge et al^ 
2008f ) JHanasoge and Callv ( 2009 ). It is now well established that slender mag- 



netic flux tubes permit the propagation of the two basic types of magnetohydro- 
dynamic waves: The longitudinal tube waves (sausage modes) with azimuthal 
wave number m = 0, which are axisymmetric, excited by pressure fluctua- 
tions. The transversal tube waves (kink modes) with m ==b 1, which are non- 
axisymmetric and describe the incompressible undulations of the tube, supported 
by magnetic tension. However, in these various papers, simplifying assumptions 
were used, such as the thin- flux-tube approxim ation, isotherma l or un stratified 
atmospheres, and the neglect of higher- m modes. ISanasoge et all ( 2008h included 
a stratified atmosphere and used the thin-flux-tube approximation to evaluate 
the scattering matrix associated with an /-mode wave interacting with a flux 
tube. The thin- flux-tube approximation works when the radius of the tube is 
smaller than it scale height which was 260 km. For this reason the largest tube 
that they considered was 160 km. Numerical treatment of flux tubes of arbitrary 
size is straightforward and robust. The appro ach is already yielding interest- 
ing re s ults in the conte x t of sunspots, see e.g. Cameron. Gizon, and Daiffallahl 

(|2007h . ll^asogd(l2008h ; lKhomenko. Collados. and Fellipd (|2008l ) J Cameron. Gizon. and Duvallfr 
(|2008h . lKhomenko et all (|2009h . 

In this paper we use numerical simulations to explore the response of vertical, 
small scale flux-tubes with radii between 200 km and 3 Mm. The paper is 
structured as follows: In Section 2 we present the background model and the 
equations describing the wave propagation. In Section 3 we present the results 
of a series of convergence tests. These are essential since we are considering 
tubes with small radii. In Section 4 we study the scattered wave field, Section 5 
is devoted to mode conversion, before discussing the results in Section 6. 



2. Simulation Model 



We use a slig htly modified version of the SLiM (Semi spectral Linear MHD) code 
described in Ca meron. Gizon. and Daiffallahl (120071 ) to propagate linear waves 



through an enhanced polytropic atmospher e. The enhanced p o ly trop ic atmo- 
sphere we use was constructed by following ICallv and Bogdanl (|l997l ) . Specifi- 



cally we follow those authors in considering a 4820 G purely vertical flux tube 
along the axis of which the pressure, density, and temperature vary as a poly- 
trope. The polytropic index is set as 1.5, the density at a depth of 1 Mm is 
set to 10 -5 g cm -3 . The sound speed is set to be equal to the Alfven speed at 
depth of 400 km. The external atmosphere is the enhanced polytrope obtained 
by requiring that the total (gas and magnetic) pressure be horizontal uniform, 
and by that the density is that of the original polytrope. Since the magnetic 
pressure is independent of height this means that away from the flux tube the 
pressure is increased by a comstant which is independent of height and thus does 
not affect the vertical hydrostatic force balance. We emphasise that this is the 
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same formulation as was used by ICallv and Bogdanl (I1997T ) and is a reasonable 



approximation for qualitatively studying wave propagation in the upper part of 
the convection zone. Figure [1] show a comparison between the vertical profiles 
of the sound speed, pressure, and acoustic velocity for the standard solar model 
S by Christensen-Dalsgaard and those on the flux-tube axis. We note that with 
our boundary conditions waves are strongly reflected from the upper boundary 
w hich restricts us to c onsid ering waves below the acoustic cut-off frequency. As 
in ICallv and Bogdanl (|l997l h we restrict ourselves to purely vertical flux tubes. 
In our case this is partly motivated by the fact that we are considering only 
small tubes (not sunspots) which are expected to be almost vertical until they 
expand into the overlying atmosphere. We thus are considering relatively small 
flux tubes of different sizes and are including the full three-dimensional geometry. 

Our code solves the linearized, ideal MHD equations in three dimensions. We 
write the magnetic, velocity, pressure, and density perturbation in terms of the 
displacement vector [£]. The background is assumed to be time independent. 
The primed quantities in the equations are the Eulerian perturbation, whereas 
the background variables are subscripted with 0. The gravitational acceleration 
[g] is constant, and there is no background steady flow. 

The equation governing the displacement vector is: 

Po^ = -W + P' 9+^(3' xB +j xB') (1) 

The perturbed quatities: density [/?], pressure [p], magnetic field [B]. Electric 
current [j] are written as functions of [£] as follow: 



p' = -V-(po€), (2) 
p' = <${p' + £'Vpo)-£'Vpo, (3) 
B' = Vx((x B ), (4) 

j' = V x B', (5) 

A flux tube that is almost evacuated at the surface and that has a radial 
profile is superposed on the background atmosphere. The magnetic flux tubes 
have a top- hat profile give n by B(r) = Bp exp(— r 4 /rp) where ro is the tubes 



radius, and B = 4820 G (ICallv and Bogdanl . Il997h . We use a local Cartesian 



geometry defined by the horizontal coordinates x and y. The vertical coordinate 
is z, and the simulation covers the height range from 0.2 Mm to 6 Mm below 
the solar surface. Note that the convention adopted here is that z incresaes with 
depth. The horizontal extent of the domain is x G [-20,20] Mm and y G [-10,10] 
Mm. The axis of the vertical flux tube passes through the point x = -7 Mm, y 
= 0. 
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Figure 1. The vertical profiles of the pressure, density and acoustic velocity for the 
background model. Red plots are the profiles for the standard solar model S by Christensen- 
Dalsgaard (with z = corresponding to the height of 400 km in model S). Blue plots are for 
the polytrope along the flux-tube axis, and the black plots are for the enhanced polytrope in 
quiet the Sun. 



At to = 0, an /-mode wave packet with a Gaussian profile in Fourier space is 
situated at xo = -20 Mm and it propagates from left to right in the x-direction. 

Specifically the initial condition is constructed from individual /-modes each 
of which has the usual (wavenumber-dependent) exponential dependence on 
height. The phase of the each eigenmode is set so that the waves are all in 
phase at x = —20 Mm. The distance from the initial condition and the sunspot 
was chosen as the result of a compromise: On the one hand it is desirable for 
the initial wavepacket to be close to the sunspot (so that the wave attenuation 
before the wave reaches the sunspot is minimized), on the other hand there is a 
relationship between the wavepacket and cross covariance from rand om sources 
( Cameron. Gizon, and Duvalli 12008 ; Gizon, Birch, and SpruitL l2010h which only 
applies outside the "near field" of the initial condition. The frequency of each 
mode is \fgk where k independently of the details of the background quiet-Sun 
model. The vertical component of the velocity, dt£, z , a t the upper surface is our 
main diagnostic, as it is closely related to the observable Doppler velocity at disk 
centre. The amplitude of the different modes in the initial condition depends on 
their eigenfrequency. The amplitudes of the surface velocity is a Gaussian with 
a central frequency and the half width of the wave packet are 3 mHz and 1.18 
mHz respectively. 

The top boundary also corresponds to that of Callv and Bogdanl ( 1997h and 
has also been previously described in ICameron. Gizon, and Duvalli ([20081 ). It 
corresponds to a free surface in the non-magnetic regions and hence is effective 
in reflecting waves. This means that the effective acoustic cut-off of the model is 
infinite. Our initial condition has very little energy (« 0.06% above 5.3 mHz). 

We have also performed a simulation without the flux tube being present, 
that is with the enhanced polytrope being used everywhere and with Bq = 0. 
This simulation acts as a reference and allows us to construct the scattered wave 
field as the difference between the simulations with and without the flux tube. 

An example of the wave interacting with a tube with a radius of 500 km is 
shown in Figure [2j For a tube of this size, the scattered wave is almost circular, 
with a significant amount of back scattering. The SLiM code assumes that the 
box is periodic, which accounts for the wave fronts coming in from the side of 
the box after 4000 seconds. It is to be noted that the diameter of the tube is 
neither small nor large with respect to the wavelength of the incoming waves, so 
that various types of tube modes will be excited. 
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Figure 2. Wave propagation at three times (1600, 3000, and 4000 seconds after the start 
of the simulation), in the presence of a 500 km magnetic flux tube. The top panel shows the 
scattered component of the wave field near the "solar surface" , the middle panel shows the full 
wave field whilst the lower panel is a vertical cut. The location and size of the tube is indicated 
by the circle at x = -7 Mm, y = 0. The point "A" is the location at which the resolution test 
analysis was performed (Section [3]). Notice that the simulation is periodic in the horizontal 
directions, which can be seen in the scatterd wave field at later times. 



3. Resolution Test 

Our simulation box has a length of 40 Mm in the direction in which the wave 
is propagating. The radii [r Q ] of the flux tubes we are interested in vary from 
200 km to 3 Mm. Furthermore, since we have used a top-hat profile [B(r) = 
Bo exp(— r 4 /ro)] the "sidewalls" of the tube are even thinner. This raises the 
question of how many Fourier modes are required in the simulation to resolve 
the flux tubes. In the vertical direction we have found that 200 grid points are 
sufficient for all the flux tubes, probably because the flux tubes have no special 
structure in the ^-direction. The tubes radii we have considered are 200 km, 
500 km, 1 Mm, and 3 Mm. We have considered three different choices for the 
number of horizontal Fourier modes: 50x25, 100x50, 200x100. The data of the 
scattered wave field at the point "A" is shown as a function of time in the top 
four panels of Figure [3l The three different resolutions are shown using different 
colors. It is apparent that the scattered wave at "A" is resolved for ro equal 
to 1 Mm and 3 Mm for all three choices of the number of Fourier modes. The 
degree of convergence in the ro = 500 km case is reasonable, especially for the 
two higher resolutions. The case where ro = 200 km is clearly a long way from 
convergence for at least the two lower resolutions. At the lowest resolution, the 
flux tube is significantly smaller than 40 Mm/50 modes = 800 km, for this reason 
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the simulated waves apparently do not even see the tube and the scattered wave 
is absent. With 100 modes, the simulation appears more reasonable but the 
amplitude is only about 1/2 of that of the simulation with 200 modes. The 
results using 200 modes should be much closer to the converged value: The 
relative resolution is similar to the case where we used 100 modes for the 500 
km tube. The four lower panels of Figure [3] show the scattered wave along the 
line y = as a function of x at a fixed time. The convergence properties are the 
same as we found for fixed position "A" . To quantify the errors we introduce a 
measure for the difference between the simulations. Again concentrating on the 
point "A", we compare the time series of the vertical velocity from either the 
low or medium-resolution simulations [v z ] with the high resolution simulation 

E f M*)-^)] 2 a. 

We can also construct a measure of the error based upon the cuts at fixed 
time (i.e. the lower four panels of Figure [3]). The error in this case is taken to be 

E*t£(*) 2 ' [> 

Figure S] shows these measures of the error in various ways. In the left panels 
we see that the error decreases rapidly as the tube increases in radius. The panels 
on the right show that increasing the resolution by a factor of two decreases this 
measure of the error by a factor of about ten. This would suggest that the error 
associated with the highest resolution (200 modes) when applied to the smallest 
flux tube (200 km radius) should be approximately 3% (i.e. a factor of ten better 
than the simulation with 100 modes which has a 30 % error). 



4. Scattering of Tubes 

The dependence of the magnitude of the scattering on radius can be seen in 
Figure [5] where the full and the scattered wave field are shown for tubes with 
radii 200 km, 500 km, 1 Mm, and 3 Mm. A vertical cut through the spot, along 
the x-axis, is shown in Figure [6j for tubes with radii 200 km, 1 Mm, and 3Mm. 
For the two larger tubes, some mode conversion is apparent. 

Naturally the simulations give us the full vector- velocity field and hence we 
can look simultaneously at all three components of the velocity field as shown 
in Figure [71 This figure makes it clear that the scattered field consists mainly 
of different mixtures of m = and m = 1 modes for different radii. From this 
figure we can see that tubes with radius between 500 km and 2 Mm are excited 
primarily in the m = mode, the upper part of the cylinder moving inward or 
outward (sausage modes). For the ^-component of velocity, tubes with radii less 
than 500 km are excited with m = ±1 dipole oscillation whereas tubes between 
500 km and 2 Mm oscillate with a mixture of m = and m = 1 modes. It is 
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Figure 3. Comparison of results at different resolutions for various sizes flux tubes. The 
top four panels show the scattered component of the wave at a fixed point ("A" in Figure [2]) 
as a function of time. The lower four panels show the scattered wave component at different 
locations along the x-axis at fixed time (corresponding to the middle panels in Figure [2]) 



SOLA: Daiffallah_revised.tex; 2 November 2010; 10:31; p. 7 



K.Daiffallah et al. 




Figure 4. The error function for scattered wave field of vertical velocity, top panel at 7 Mm 
from tube center along y = 0, bottom panel at time = 75 minutes along y = 0. Curve with 
stars at left correspond to the difference low vs high resolutions, curve with cross is for the 
difference medium vs high resolutions. The right curves shows the variation of the error for 
the different radius when the horizontal resolution increases. 
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Figure 5. Wave propagation through tubes of different sizes. The images are from snapshots 
taken 3000 seconds after the start of the simulation. The upper four panels show vertical 
velocity scattered wave field (at the top) and the full wave field (at the bottom) for tubes with 
radii 200 km and 500 km, the lower four panels are for tubes with radii 1 Mm and 3 Mm. The 
circles indicate the sizes of these tubes. The simulation is periodic in the horizontal directions 
which can be seen in the scattered wave field. 
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Figure 6. The x-component of the velocity near the photosphere; a vertical cut passing 
through the centre of the tubes 200 km, 1 Mm, 3 Mm. All panels are from a snapshot at t = 
2000 seconds. 
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Figure 7. Scattered wave field for different tubes radii taken 3000 seconds after the start of 
the simulation. Radius from top to bottom: 200 km, 500 km, 3 Mm. Velocity component from 
left to right: V x , V y , V z . Note that different gray-scales are used for each frame so that only 
the structure of the scattered wavefield is significant 



more difficult to identify the oscillation of 3 Mm tube radius, which presumably 
has substantial energy in higher m components. 

To procede further we consider the displacement of the tube with respect to 
its surroundings. The idea is illustrated in Figure [51 Thinner flux tubes (radius 
smaller than 500 km) are effectively coupled to the motion of the surrounding 
plasma by the drag force, while larger tubes (radius larger than 500 km) can move 
relative to the surrounding plasma owing to the action of the magnetic force, sim- 
ilar t o a flexible solid body immersed in a fluid ([Solanki, Inhester, and Scfmssler , 
2006[ ). This result was expected, but is still nice to see in the simulations. 

To better understand the m modes that have been excited in the tube, we 
look at the total ^-displacement not only along axis of the tube but also near its 
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Figure 8. Horizontal displacement in arbitrary units in the x-direction of the tube axis and 
the background as a function of the depth [z]. The curves are taken 2000 seconds after the 
start of the simulation for 200 km tube radius, left, and 3 Mm tube radius, right. 
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Figure 9. Horizontal displacement in arbitrary units in the x-direction of the tube axis 
(dashed line at the center) and the tube cross sections (full line right and left) as a function of 
the depth [z]. The curves are taken 1600 seconds after the start of the simulation for different 
tubes radii: 200 km, 500 km, 1 Mm, and 3 Mm respectively. 



SOLA: Daiffallah_revised.tex; 2 November 2010; 10:31; p. 11 



K.Daiffallah et al. 




Figure 10. The x-component of the perturbation of the magnetic field introduced by the 
passage of the wave, shown in the x — z (left) and y — z (right) planes passing through the axis 
of the sunspot. The downward-propagating oscillations correspond to the slow- acoustic mode 
propagating along the fieldlines. This figure is for t = 2280 seconds and the 3 Mm flux tube. 

edges. This is shown in Figure [9J Consistent with the above description, we see 
that the perturbations are coherent across small tubes. The change in behavior 
around 500 km corresponds to the fact that this is near the characteristic length- 
scale of the / mode (g/uj 2 ) of 771 km for the 3 mHz centered wave packet. For 
larger tubes, higher-m modes begin to dominate because there is little coherence 
of the incoming wave across the tube. 



5. Mode Conversion 



As in ICallv and Bogdanl ( 1997 ). we find that the magnetic flux tubes convert 



some of the incoming waves into slow-magnetoacoustic modes which propagate 
along the field-lines. This can be seen in Figure [8j where the horizontal dis- 
placement of the flux with respect to the background is shown. The vertical 
oscillations shown here are one signature of downward propagating slow mode 
waves. The magnetic nature of these waves can be easily seen in Figure [TUl where 
a component of the magnetic-field perturbation is shown. We comment that the 
mode conversion removes energy from the incoming waves. 



6. Conclusions 

We have carried out a series of numerical simulations of plane- wave propagation 
through flux tubes of different sizes. In the first part of this study, we tested the 
convergence in order to determine the minimum-sized tube that we could reliably 
simulate with the resolutions that we had available. We then investigated how 
different-sized tubes interacted with an incoming /-mode wavepacket. Different 
scattered wave-field patterns were observed for flux tubes with different radii. 
In agreement with previous suggestions, we found that when the flux tube was 
small compared with the wavelength of the incoming wave, mainly the m = 1 
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kink modes were excited. For mid-ranged tubes both the m = 1 kink and m = 
sausage modes were excited, and for large tubes numerous m and radial order 
modes were excited. Our results demonstrate that numerical simulations are an 
efficient, robust, and straightforward method to treat the action of flux tubes on 
waves regardless of their size. The simulations are a rich source of information, 
which can be used to better understand the observations in future studies. 
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